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Abstract.  Theoretical  understanding  and  numerical  modeling  of  atmospheric  moist 
convection  still  pose  great  challenges  to  meteorological  research.  A  Direct  Numerical  Sim¬ 
ulation  of  a  single  cumulus  cloud  is  beyond  the  capacity  of  today’s  computing  power.  The 
use  of  a  Large  Eddy  Simulation  in  combination  with  semi-implicit  time-integration  and 
adaptive  techniques  offers  a  significant  reduction  of  complexity. 

This  paper  presents  a  first  step  towards  an  efficient  simulation  of  a  single  cloud.  So  far 
this  work  is  restricted  to  dry  flow  in  two-dimensional  geometry  without  subgrid-scale  mod¬ 
eling.  The  compressible  Euler  equations  are  discretized  using  a  discontinuous  Galerkin 
method  introduced  by  Giraldo  and  Warburton  in  2008.  Time  integration  is  done  by  a 
semi-implicit  backward  difference  formula.  This  paper  represents  the  first  application  of  a 
triangular  discontinuous  Galerkin  method  with  h-adaptive  grid  refinement  for  nonhydro¬ 
static  atmospheric  flow.  This  refinement  of  our  triangular  grid  is  implemented  with  the 
function  library  AMATOS  and  uses  an  efficient  space  filling  curve  approach. 

Validation  through  different  test  cases  shows  very  good  agreement  between  the  current 
results  and  those  from  the  literature.  For  comparing  different  adaptivity  setups  we  devel¬ 
oped  a  new  qualitative  error  measure  for  the  simulation  of  warm  air  bubbles.  With  the 
help  of  this  criterion  we  show  that  the  simulation  of  a  rising  warm  air  bubble  on  a  locally 
refined  grid  can  be  four  times  faster  than  a  similar  computation  on  a  uniform  mesh  while 
still  producing  the  same  accuracy.  Remarkably  only  5%  of  the  total  CPU  time  is  used  for 
adapting  the  grid  after  each  time-step. 
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1  INTRODUCTION 

Significant  progress  in  numerous  areas  of  scientific  computing  comes  from  the  steadily 
increasing  capacity  of  computers  and  the  advances  in  numerical  methods.  An  example  is 
the  simulation  of  the  Earth’s  atmosphere,  which  has  proven  to  be  extremely  challenging 
owing  to  its  multiscale  and  multi-process  nature.  Even  with  todays  computers  it  is  impos¬ 
sible  to  explicitly  represent  all  scales  and  all  processes  involved.  To  overcome  this  difficulty 
one  resorts  to  empirically-based  closure  approaches  —  called  “parameterisations”  —  that 
try  to  capture  the  unresolved  aspects  of  the  problem.  Needless  to  say,  this  introduces 
errors. 

An  application  with  high  practical  relevance  is  numerical  weather  prediction.  Generally, 
its  skill  has  improved  considerably  over  the  past  few  decades,  and  a  significant  portion 
of  this  improvement  has  been  attributed  to  the  increased  computing  power  and  refined 
numerical  methods1.  A  notable  exception  to  this  general  development  is  the  forecast 
of  precipitation,  where  the  progress  has  been  almost  non-existent2.  The  reason  for  this 
state  of  affairs  is  likely  due  to  the  fact  that  most  processes  leading  to  precipitation  are 
parameterized  rather  than  explicitly  simulated  in  today’s  weather  prediction  models.  In 
this  paper  we  shall  develop  and  present  a  new  numerical  model  that  is  specifically  tailored 
to  investigate  one  of  these  processes  in  detail. 

The  paper  is  organized  as  follows.  First,  in  section  2  we  motivate  the  new  model 
development  by  explaining  the  meteorological  problem.  In  section  3  we  then  present  the 
numerical  methods  used  in  our  work.  This  includes  the  discontinuous  Galerkin  method 
for  the  spatial  discretization,  a  semi- implicit  method  for  the  time  integration,  and  a  space 
hlling  curve  approach  for  the  adaptive  grid  management.  In  section  4  we  apply  our  code 
to  three  test  cases  from  the  literature.  Section  5  provides  some  tests  concerning  the 
accuracy  of  the  adaptive  mesh  refinement.  The  paper  ends  with  a  summary  and  outlook 
in  section  6. 

2  METEOROLOGICAL  MOTIVATION 

A  single  cumulus  cloud  can  be  considered  as  a  prototypical  basic  element  of  atmo¬ 
spheric  moist  convection.  The  cloud  rises  through  the  environmental  air  owing  to  its 
positive  buoyancy  (fig.  1).  Upward  motion  of  the  cloud  (thick  blue  arrow)  is  associated 
with  downward  motion  in  a  thin  shell  surrounding  the  rising  cloud  (thin  blue  arrows)3. 
This  induces  wind  shear  at  the  cloud-environment  interface,  leading  to  Kelvin-Helmholtz 
instabilities  which  eventually  result  in  turbulence.  The  ensuing  mixing  between  moist 
cloudy  and  dry  environmental  air  leads  to  evaporation  of  cloud  droplets.  This  cools  the 
parcel  resulting  in  negative  buoyancy  corresponding  to  a  downward  force  (red  arrows). 
This  process  is  aptly  called  “buoyancy  reversal”  5,6,7 . 

Early  indications  for  the  significance  of  buoyancy  reversal  for  cloud  dynamics  stem  from 
the  laboratory  experiments  of  Johari8.  He  introduced  buoyancy  reversal  in  his  watertank 
experiments  with  the  help  of  chemical  reactions  occurring  in  the  mixing  region  between 
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Figure  1:  Illustration  of  buoyancy  reversal.  The  blue  arrows  demonstrate  the  mean  flow 
of  a  rising  cloud,  the  black  arrows  represent  turbulence  produced  by  Kelvin-Helmholtz 
instability  and  the  red  arrows  illustrate  buoyancy  reversal.  For  further  explanation  we 
refer  to  the  text. 


two  fluids.  Johari  found  that,  depending  on  the  strength  of  the  buoyancy  reversal,  the 
morphology  of  the  cloud  development  could  be  vastly  different. 

Similar  results  were  found  in  highly  idealized  numerical  two-fluid  experiments  by  Gra- 
bowski1  in  1995.  These  simulations  started  with  two  fluid  layers,  one  on  top  of  the 
other.  Convergence  was  imposed  in  the  lower  layer  thus  leading  to  a  rising  plume.  The 
simulation  was  done  twice:  with  and  without  buoyancy  reversal.  The  buoyancy  reversal 
was  implemented  as  an  additional  downward  force  at  the  interface  between  the  two  fluids. 
Initially  the  two  simulations  looked  very  similar  but,  after  further  development,  increasing 
differences  could  be  seen. 

These  preliminary  investigations  suggest  that  buoyancy  reversal  has  an  important  im¬ 
pact  on  cloud  dynamics  and,  hence,  on  the  formation  of  precipitation.  However,  owing  to 
their  idealized  nature  it  is  not  possible  to  draw  any  firm  conclusion  about  real  clouds.  On 
the  other  hand,  numerical  weather  prediction  models  and  even  so-called  “cloud  resolving 
models”  are  not  able  to  explicity  simulate  the  processes  relevant  for  buoyancy  reversal 
due  to  their  coarse  spatial  resolution10.  It  is  here  that  we  want  to  take  a  step  forward  by 
developing  a  new  numerical  model  that  is  specifically  designed  to  deal  with  the  mixing 
processes  at  the  cloud  boundary. 

A  Direct  Numerical  Simulation  (DNS)  would  require  a  resolution  of  about  1  mm  in 
each  direction  in  order  to  properly  resolve  all  dynamical  scales10.  In  three  dimensions 
this  amounts  to  some  1024  grid  points,  which  is  beyond  the  capacity  of  today’s  computing 
power.  We,  therefore,  resort  to  Large  Eddy  Simulation  (LES)  in  combination  with  a  semi- 
implicit  time  integration  method.  The  use  of  an  adaptive  technique  will  offer  a  significant 
reduction  in  numerical  expense,  as  it  allows  us  to  focus  attention  to  the  cloud-environment 
interface,  which  is  the  region  where  mixing  and  buoyancy  reversal  takes  place. 
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3  NUMERICAL  METHODS 

In  this  section  we  present  the  numerical  methods  that  are  used  in  our  work.  The  choice 
of  the  numerical  method  requires  some  a  priori  knowledge  of  the  application  for  which  the 
numerical  model  is  designed.  Measurements  in  real  clouds  have  shown  a  large  variety  of 
behavior9:  there  are  clouds  with  numerous  steep  gradients  in  the  interior  (see  for  example 
the  liquid  water  content  in  figure  15  of  that  reference9).  On  the  other  hand,  smaller  clouds 
often  have  a  fairly  smooth  interior  with  discontinuities  mostly  at  the  boundary  of  the  cloud 
(fig.  13  of  Damiani  et  al.9).  It  is  the  latter  which  we  intend  to  simulate.  For  such  clouds 
a  discontinuous  Galerkin  (DG)  discretization  in  combination  with  a  semi-implicit  time- 
integrator  should  be  the  best  choice;  the  reason  for  this  choice  is  due  to  the  high-order 
accuracy  and  robustness  in  handling  discontinuities  of  the  DG  method  as  well  as  the  large 
time-steps  allowed  by  the  semi-implicit  method.  For  avoiding  artificial  oscillations  an 
artificial  viscosity  is  used.  Before  describing  these  methods  in  the  following  subsections 
we  first  introduce  the  equations  that  we  use  in  our  simulations. 

In  the  current  work  the  fully  compressible  Euler  equations  are  used.  For  the  dry  case 
we  use  the  following  set  of  equations  (see  equation  set  2  in  Giraldo  and  Restelli13): 


!+v.(pU) 

=  0, 

(i) 

+  V  •  (pu  ®  u  +pl2) 

=  ~pg  k, 

(2) 

dp9 

at  +v'(p9u) 

=  0, 

(3) 

T  T 

where  the  variables  are  (p,  pu,  p9 )  ,  p  is  the  density,  u  =  (u,  w)  is  the  velocity  held  and 
9  is  the  potential  temperature.  Furthermore  we  denote  the  gravitational  constant  with  g, 
the  tensor  product  by  the  identity  matrix  in  M2  by  1 2 ,  the  unit  vector  in  the  vertical 
direction  with  k. 

Pressure  p  in  eq.  (2)  is  given  by  the  equation  of  state: 


P  =  Po 


(4) 


with  a  constant  reference  pressure  po,  the  gas  constant  R  =  cv  —  cv  and  the  specific  heats 
for  constant  pressure  and  volume,  cp  and  cv.  Potential  temperature  9  is  defined  by 


9 


(5) 


with  temperature  T.  This  definition  can  be  illustrated  in  the  following  way:  if  we  consider 
dry  air  with  temperature  T  and  pressure  p,  potential  temperature  9  is  given  by  the 
temperature  which  the  air  would  have  when  being  transported  adiabatically  to  a  place 
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with  pressure  p0.  In  our  model  we  use  potential  temperature  as  a  variable  because  this 
simplifies  the  extension  to  moist  air  in  future  research. 

Atmospheric  flow  is  often  approximately  in  hydrostatic  balance,  which  is  defined  by 


dp 

dz 


=  - pg ■ 


(6) 


This  balance  can  produce  numerical  instabilities  because  the  remaining  terms  in  the  ver¬ 
tical  component  of  eq.  2  are  much  smaller  than  the  two  terms  of  the  hydrostatic  balance 
(6).  To  avoid  this  instability  we  introduce  the  mean  states  p,  p  and  9  that  are  in  hydro¬ 
static  balance.  The  mean  state  of  pressure  p  is  defined  by  p  —  p(p,6).  The  deviation  of 
the  variables  from  the  mean  state  is  denoted  by  p'  =  p  —  p,  9' =  6  —  0  and  p'  =  p  —  p.  By 
this  procedure  the  set  of  equations  (1)  -  (3)  can  be  written  as 


dpu 

~dT 


dp' 

V  ■  {pu®u-\-p'I2)  =  -  p'  g  k, 

dp6' 


dt 


+  ~V  ■  (p  9  u)  =  0. 


To  discretize  these  equations  in  space  we  introduce  the  commonly  used  notation 

f  +  V-F  (,)  =  S(,). 

with  the  vector  q  =  (p,  pu,  p0)T ,  the  source  function 

0 

s(q)  =  |  -p'fi'k  |  , 

0 


(7) 

(8) 

(9) 

(10) 


(11) 


and  the  flux  tensor 


pu 

F  (q)  —  (  pu®u  +  p'  I2 
p6u 


(12) 


Equation  (10)  is  discretized  using  the  discontinuous  Galerkin  method  which  we  now  de¬ 
scribe  in  detail. 


3.1  Discontinuous  Galerkin  Method 

In  our  work  we  use  a  discontinuous  Galerkin  method  based  on  the  strong  formulation 
using  the  Rusanov  flux  at  the  cell  interfaces.  Furthermore,  we  consider  a  two  dimensional 
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triangular  mesh;  the  extension  to  a  full  three  dimensional  method  will  remain  a  task  for 
the  future  but  we  envision  using  tetrahedra  for  this  task.  The  triangular  discontinuous 
Galerkin  method  used  in  our  work  is  described  by  Giraldo  and  Warburton11  for  the  case 
of  the  shallow  water  equations.  Despite  a  different  definition  of  conserved  variables  q , 
flux  tensor  F  (q)  and  source  function  S  (q),  eq.  (10)  remains  unchanged.  Therefore,  we 
repeat  in  this  paper  only  the  main  ideas  of  the  discretization. 

We  start  by  multiplying  eq.  (10)  with  a  test  function  ip,  integrating  over  an  arbitrary 
element  Vte  and  bringing  the  spatial  derivative  in  front  of  the  test  function  with  integration 
by  parts.  Replacing  the  flux  in  the  boundary  terms  by  a  numerical  flux  F*  leads  to  the 
following  equation  for  the  numerical  solution  qN : 


dq 


N 


dt 


—  Ftv-V  —  Sn)  ^  (x)  dx 


ip  (x)  n  ■  F^dx, 


(13) 


where  Te  is  the  boundary  of  element  fie,  n  is  the  outward  pointing  unit  normal  vector  on 
Te,  F tv  =  F  ( qN )  and  Sjv  =  S  ( qN ).  Applying  again  integration  by  parts  gives  the  strong 
formulation 12 


+  V  •  Fjv 


if  ( x )  dx 


tp  (*)  n  •  ( Fn  -  F*n)  dx. 


Now  we  introduce  an  expansion  by  polynomial  basis  functions 


(14) 


mn 

Qn  (*)  =  ^2ipj(x)qj 

3= 1 


(15) 


and  assume  that  the  test  function  ip  can  be  written  as  a  linear  combination  of  the  basis 
functions.  Using  Einstein’s  sum  convention  we  get 


dqj 

dt 


3pi  (V  •  Fn  -  SN)  dx  +  /  xpih  dx  •  (F 


N 


N) 


(16) 


where  3pi  =  1  -0fc  with  the  mass  matrix  M. ^  =  fn  ^^dx.  For  the  sake  of  simplicity,  we 

did  not  write  the  dependence  on  x  of  the  basis  functions  although  it  should  be  understood 
that  the  basis  functions  depend  on  the  spatial  coordinates. 

The  integrals  in  eq.  (16)  are  evaluated  using  high  order  quadrature  rules11.  F*  is 
discretized  as  the  Rusanov  flux,  given  by: 

F«G[F(gt)+FM)- An  (17) 

with  the  maximum  wave  speed  A  =  max  (|u|  +  a,  |u|  —  a)  where  a  is  the  speed  of  sound. 
If  the  normal  vector  n  of  element  h2e  is  pointing  to  the  right,  q ^  is  the  left  limiting  value 
of  qN  and  q^f  is  the  right  limiting  value. 
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So  far  we  have  derived  a  discontinuous  Galerkin  discretization  for  our  set  of  equations 
(10).  At  this  point,  the  right  hand  side  of  eq.  (16)  is  known  and  we  can  integrate 
the  equation  in  time.  This  can  be  done  either  by  an  explicit  or  implicit  method.  For  an 
explicit  method  we  implement  a  third  order  Runge-Kutta  method  of  Cockburn  and  Shu15. 
Because  of  the  fast  sound  and  gravity  waves  this  explicit  time-integration  is  restricted  to 
a  very  short  time-step.  As  explained  before  we  are  not  interested  in  simulating  these  fast 
waves  accurately;  therefore,  we  also  use  a  semi-implicit  time-integrator  as  presented  in 
the  next  subsection. 


3.2  Semi-Implicit  Time  Integration 

The  semi- implicit  time  integration  is  implemented  in  a  similar  fashion  to  the  approach 
of  Restelli  and  Giraldo16,17.  The  main  difference  is  that  we  use  potential  temperature 
instead  of  total  energy  as  the  fourth  variable. 

The  full  nonlinear  operator  J\f  ( q )  is  given  in  our  notation  by 


AT(q)  =  -V-F  (q)  +  S(q). 
For  the  semi-implicit  approach  we  define  an  operator  C  by 


Cq  = 


/  V-(p«)  \ 

dp' /dx 
dp'  jdz  +  g  p' 

\V-(Spu)  J 


(18) 


(19) 


where  ft  is  given  by  ft  =  p9' .  The  operator  C  is  linear  in  the  conserved  variables 

(p,  pu,  p9)  and  p'  is  a  linearized  version  of  pt  with  respect  to  the  conserved  variables.  As 
explained  by  Restelli19  the  operator  C  is  responsible  for  the  fast  moving  sound  and  gravity 
waves  and,  therefore,  must  be  integrated  implicitly.  This  splitting  is  done  by  writing 

=  {A f  {q)  -  Cq}  +  Cq.  (20) 

For  discretizing  (20)  in  time,  we  use  a  backward  difference  formula  of  order  2,  that  leads 
to 


1  1  1 

X  =  A  W  (<?”)  -  £«”]  +  ■c9“+1  (21) 

^  m=— 1  m= 0 

with  a_i  =  1,  cko  =  4/3,  aq  =  —1/3,  7  =  2/3,  (3q  =  2  and  =  —1.  We  rewrite  this 
equation  collecting  all  terms  with  qn+1  and  get 


1 

[1  -  7  At  C]  qn+1  =  qex~l^tJ2  PmC  qn-m,  (22) 

m= 0 
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where 


1  1 

<T  =  amqn-m  +  7  At  ]T  fimN  ( qn~m )  (23) 

m= 0  m= 0 

is  an  explicit  predictor  that  has  to  be  calculated  first.  Solving  the  linear  system  of  equa¬ 
tions  (22)  (e.g.,  with  GMRES)  gives  the  implicit  corrector.  For  details  on  this  solution 
strategy  the  reader  is  referred  to1'  and18. 

3.3  Mesh  Refinement  with  Space  Filling  Curve  Approach 

As  explained  in  section  2  we  expect  steep  gradients  at  the  boundary  of  the  cloud.  For 
increasing  the  numerical  resolution  in  these  regions  we  use  h-adaptive  mesh  refinement. 
This  is  managed  with  the  function  library  AMATOS23.  The  main  advantage  of  this 
function  library  is  that  it  handles  the  entire  h-adaptive  mesh  refinement.  Furthermore  it 
orders  the  unknowns  very  efficiently  by  using  a  so-called  space  filling  curve  approach.  For 
further  information  we  refer  to  the  publication  of  Behrens  et  al.23. 

The  only  modification  that  was  necessary  for  our  work  was  the  calculation  of  the  new 
values  at  the  grid  points  when  elements  are  refined  or  coarsened.  This  is  quite  straight 
forward.  We  simply  evaluate  the  old  polynomials  at  the  positions  of  the  new  degrees  of 
freedom. 

The  refinement  criterion  used  for  the  results  presented  in  this  paper  is  given  by: 

\0\x,t)\  >  max  (\9'(x,  t)|)  / 10.  (24) 

X 

Wherever  this  condition  is  fulfilled  the  mesh  is  refined  until  it  reaches  a  specified  finest 
resolution.  In  the  rest  of  the  domain  the  grid  is  coarsened  until  it  reaches  a  specified 
coarsest  resolution  without  modifying  the  resolution  in  the  refinement  region.  The  tran¬ 
sition  between  fine  and  coarse  meshes  is  given  by  the  conformity  of  the  grid.  To  avoid 
small  scale  structures  moving  into  a  region  with  a  coarse  mesh  we  add  two  to  four  rows  of 
fine  elements  to  the  refinement  region.  For  more  realistic  simulations,  we  intend  to  utilize 
more  involved  refinement  criteria,  which  consider  physical  behavior,  gradients,  and  other 
indicators. 

3.4  Artificial  viscosity 

For  avoiding  artificial  oscillations  we  use  a  so  called  artificial  viscosity.  This  is  done  by 
adding  a  diffusion  term  V  ■  (/tpVu)  to  the  right  hand  side  of  eq.  (2)  and  V  ■  (pp’VO') 
to  the  right  hand  side  of  eq.  (3)  with  the  viscosity  parameter  p.  These  terms  represent 
the  simplest  variational  multiscale  approach  one  could  use. 

For  discretizing  the  second  order  derivates  produced  by  the  artificial  viscosity  we  use 
a  local  discontinuous  Galerkin  method14.  That  means  that  the  order  of  the  derivatives  is 
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reduced  by  introducing  the  following  new  variables 

a  =Vm,  (25) 

P  =V0.  (26) 

Similar  to  eq.  (16)  we  get 

cti  =  —  /  Uj'V'ifjjdx  +  /  '(/y^fidx  (ttj  —  u*)  ,  (27) 

Pi  =  -  [  i>i  dj'Vi/jjdx.  +  f  't^i'ipjhdx  (dj  ~  0*)  .  (28) 

J  Qe  ^  Te 

As  these  viscosity  terms  do  not  describe  a  flow  in  a  certain  direction  (as  in  the  case  of 
the  advection  terms)  we  use  the  following  numerical  flux  for  the  viscosity  terms  in  F*N  as 
well  as  for  u*  and  9*: 

Ft,c(<m  =  )(F(«J)+F  («£;))■  (29) 

Note  that  this  is  not  the  only  possibility  for  discretizing  the  second  order  operators;  for 
other  choices  see  Shahbazi  et  al.20. 

For  the  density  current  test  case  of  Straka  et  al.22  we  use  the  viscosity  parameter  of 
/j,  =  75m2/s  as  defined  by  Straka  et  al.22.  The  other  results  presented  in  this  paper  are 
obtained  with  a  viscosity  parameter  of  /i  =  0.1m2/s.  We  found  that  this  artificial  viscosity 
is  very  well  suited  for  avoiding  artificial  oscillations  for  the  numerical  resolutions  used  in 
this  paper. 

4  VALIDATION 

For  the  validation  of  our  new  numerical  model  we  considered  three  test  cases  that  are 
relevant  for  atmospheric  convection.  These  test  cases  are  a  small  cold  air  bubble  on  top 
of  a  large  warm  air  bubble  from  Robert21,  a  density  current  from  Straka  et  al.22,  and  a 
smooth  warm  air  bubble  from  Giraldo  and  Restelli13.  In  these  cases  no  exact  solution 
exists  but  we  can  compare  our  results  with  those  from  the  literature. 

4.1  Small  Cold  Air  Bubble  on  Top  of  Large  Warm  Air  Bubble 

The  first  test  case  that  we  consider  is  a  small  cold  air  bubble  on  top  of  a  large  warm  air 
bubble  in  a  domain  of  1km x  lkm.  This  test  case  was  introduced  by  Robert21  in  1993.  The 
background  state  has  a  constant  potential  temperature  of  9  —  300  K.  Both  bubbles  have 
a  Gaussian  profile  in  9’ .  The  warm  air  bubble  has  an  amplitude  of  0.5  K,  the  amplitude 
of  the  cold  air  bubble  is  0.17  K.  The  initial  conditions  are  chosen  identically  to  those  in 
the  publication  of  Robert21. 

In  our  test,  however,  we  introduce  two  modifications.  First,  we  use  a  slightly  different 
resolution;  the  shortest  element  edge  in  our  simulation  has  a  length  of  11m.  In  combination 
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with  third  order  polynomials  the  number  of  degrees  of  freedom  corresponds  to  a  first  order 
method  with  a  resolution  of  about  4m.  We  call  this  reduced  value  “effective  resolution”. 
With  4m  it  is  slightly  smaller  than  5m  of  Robert21. 

Figure  2  shows  our  result  for  this  test  case.  The  mesh  is  continuously  adapted  to  the 
position  of  the  temperature  anomaly.  For  this  reason  the  fine  mesh  follows  the  bubbles 
very  nicely.  By  comparing  our  result  with  the  corresponding  figure  of  Robert21  one  can 
see  that  the  results  agree  very  well.  After  600s  the  position  and  shape  of  the  warm  air  is 
still  almost  identical  to  the  corresponding  plot  of  Robert21.  Even  the  smaller  vortices  on 
the  right  hand  side  of  the  domain  are  very  similar  to  the  result  of  Robert. 

4.2  Density  Current 

A  second  test  case  is  a  density  current  initialized  by  a  cold  air  bubble  with  a  cosine 
profile  and  an  amplitude  of  15  K  in  9'  (figure  3).  This  test  case  was  introduced  by  Straka 
et  al.22.  The  viscosity  of  fi  =  75m2/s  is  identical  to  the  setup  of  Straka  et  al.22.  As  in 
the  previous  test  case  we  use  third  order  polynomials.  The  shortest  element  edge  has 
a  length  of  50m  in  our  computation.  This  leads  again  to  an  effective  resolution  that  is 
slightly  smaller  than  the  smallest  resolution  of  Straka  et  al.  with  25m.  Again  we  see  no 
differences  between  our  result  and  the  results  in  the  literature.  The  position  of  the  density 
current  and  the  shape  of  the  Kelvin-Helmholtz  rotors  agrees  very  well  with  the  result  of 
Straka  et  al.22  throughout  the  whole  simulation. 

4.3  Smooth  Warm  Air  Bubble 

As  a  third  test  case  we  computed  the  rising  thermal  bubble  introduced  by  Giraldo  and 
Restelli13  (test  case  2)  in  2008.  ft  is  a  single  warm  air  bubble  with  a  cosine  profile  in  9' . 
As  in  the  test  case  of  Robert21  the  domain  has  an  extent  of  1km  in  each  direction  and  the 
bubble  has  an  amplitude  of  0.5  K.  We  use  the  same  parameters  as  in  the  publication  of 
Giraldo  and  Restelli13  except  an  additional  artificial  viscosity  and  a  modified  numerical 
resolution.  In  our  case  the  shortest  element  edge  has  a  length  of  11m.  By  using  third 
order  polynomials  our  effective  resolution  is  about  4m.  This  is  significantly  larger  than  all 
the  effective  resolutions  used  by  Giraldo  and  Restelli.  As  they  use  10th  order  polynomials 
one  has  to  divide  their  resolution  by  a  factor  of  three  to  get  a  comparable  number  of 
degrees  of  freedom  as  in  our  simulation.  This  explains  that  they  get  almost  no  artificial 
oscillations  even  without  using  artificial  viscosity. 

As  in  the  previous  test  cases  there  are  no  obvious  differences  between  our  results  and 
those  from  the  literature.  The  position  and  shape  of  the  bubble  seems  to  be  identical  to 
the  result  of  Giraldo  and  Restelli13.  This  gives  us  confidence  that  our  code  is  error-free. 
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Figure  2:  Small  cold  air  bubble  on  top  of  a  large  warm  air  bubble  as  introduced  by 
Robert21.  The  contour  lines  show  the  deviation  of  the  potential  temperature  from  the 
background  state  and  the  gray  lines  show  the  adaptively  refined  triangular  mesh  used  in 
onr  simulation.  The  contour  values  are  from  -0.05  K  to  0.45  K  with  an  interval  of  0.05  K. 
For  avoiding  artificial  oscillations  we  use  a  constant  artificial  viscosity  of  /i  —  0.1m2/s. 
For  the  time-integration  we  used  the  explicit  Runge-Kutta  method.  Semi-implicit  time- 
integration  produces  the  same  results. 


11 


Andreas  Miiller,  Jorn  Behrens,  Francis  X.  Giraldo  and  Volkmar  Wirth 


x/km 


Figure  3:  Density  current  as  introduced  by  Straka  et  al.22.  As  in  figure  2  the  contour  lines 
show  the  deviation  9'  of  the  potential  temperature  6  from  the  background  state  6  and 
gray  lines  show  the  adaptively  refined  triangular  mesh.  Contour  values  are  from  0.5  K  to 
14.5 K  with  an  interval  of  IK.  Because  of  the  fairly  strong  viscosity  of  //  =  75m2/s  no 
values  outside  this  ranee  exist. 
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Figure  4:  Rising  thermal  bubble  introduced  by  Giraldo  and  Restelli13.  As  in  figure  2  the 
contour  lines  show  the  deviation  9'  of  the  potential  temperature  6  from  the  background 
state  6  and  gray  lines  show  the  adaptively  refined  triangular  mesh.  Contour  values  are 
from  0.025  K  to  0.45  K  with  an  interval  of  0.025  K.  Because  of  the  artificial  viscosity  of 
H  =  0.1m2/s  there  exist  no  negative  values  and  no  values  larger  than  0.45  K. 
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5  SENSITIVITY  STUDIES 

One  important  question  for  each  adaptive  numerical  model  is:  how  accurate  is  the 
adaptive  method?  For  comparing  different  adaptivity  setups  we  introduce  a  new  error 
measure  for  the  simulation  of  warm  air  bubbles.  In  this  section  we  start  with  describing 
this  criterion  that  is  used  later  for  some  sensitivity  studies.  These  studies  include  a 
comparison  between  a  simulation  on  an  adaptive  mesh  with  a  simulation  on  a  uniform 
mesh  and  a  sensitivity  study  concerning  the  size  of  the  refinement  region. 

5.1  Comparison  Criterion 

For  comparing  different  adaptivity  setups  we  had  to  find  a  criterion  that  gives  at  least 
an  approximation  of  the  accuracy  of  an  adaptive  simulation.  As  explained  before  no  exact 
solution  for  the  warm  air  bubble  test  case  is  known.  Hence  it  is  impossible  to  calculate 
the  error  of  those  simulations.  Furthermore  the  results  for  increasing  numerical  resolution 
do  not  converge  as  more  and  more  turbulence  is  resolved.  Instead  we  sought  a  qualitative 
measure  for  the  numerical  errors.  For  this  purpose  we  use  the  rising  thermal  bubble 
of  Giraldo  and  Restelli13  as  in  section  4.3.  But  this  time  we  continue  the  simulation 
much  longer  until  numerical  errors  become  clearly  visible.  After  900  seconds  we  get  the 
results  shown  in  figure  5  for  three  different  resolutions  with  a  constant  artificial  viscosity 
of  /i  =  0.1m2/s.  This  figure  suggests,  that  the  results  converge.  Furthermore  one  can  see 
that  vortices  which  arise  due  to  the  Kelvin- Helmholtz  instabilities  at  the  left  bottom  part 
of  the  bubble  strongly  depend  on  numerical  resolution.  This  is  not  surprising.  A  similar 
situation  can  be  seen  in  the  case  of  shear-instability  of  a  horizontal  shear-flow.  The  exact 
solution  of  a  perfectly  horizontal  shear-flow  is  just  stationary24.  A  perfectly  horizontal 
shear-flow  will  remain  perfectly  horizontal  as  long  as  there  is  no  perturbation.  Some  kind 
of  perturbation  is  necessary  for  making  the  instability  of  the  flow  visible.  In  a  similar  way 
we  expect  that  a  perfectly  circular  warm  air  bubble  in  an  infinitely  large  domain  should 
never  develop  Kelvin-Helmholtz  instabilities.  But  even  a  tiny  perturbation  will  grow  in 
time  and  becomes  visible  for  a  long  enough  simulation  time.  In  our  numerical  simulations 
the  perturbation  is  given  only  by  numerical  errors.  So  we  can  compare  the  accuracy  of 
different  simulations  by  comparing  the  time  when  the  instability  at  the  bottom  left  side 
of  the  bubble  develops. 

This  criterion  is  not  suitable  for  comparing  different  numerical  models.  The  onset  of 
Kelvin-Helmholtz  instability  depends  also  on  diffusion.  The  larger  the  diffusion  the  later 
the  instability  can  be  seen.  A  numerical  model  with  strong  numerical  diffusion  shows 
less  instability  but  is  not  more  accurate.  Therefore  the  artificial  viscosity  has  to  be  large 
compared  to  possible  differences  in  the  numerical  diffusion  of  the  simulations  that  are 
compared. 

A  further  remark  concerning  figure  5.  At  the  top  left  corner  of  the  domain,  one  can  see 
a  vortex  that  is  almost  independent  of  the  numerical  resolution.  This  vortex  is  caused  by 
Kelvin-Helmholtz  instability  too,  but  in  this  case  the  instability  is  initiated  by  the  solid 
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Figure  5:  Rising  thermal  bubble  as  in  figure  4  for  three  different  resolutions  with  a 
constant  artificial  viscosity  of  n  —  0.1m2/s.  Shown  is  only  the  left  half  of  the  bubble. 
Color  shading  indicates  the  potential  temperature  deviation  from  the  background  state. 
The  adaptively  refined  triangular  mesh  is  shown  by  the  black  lines.  The  mesh  is  not  shown 
where  the  refinement  criterion  (24)  is  fulfilled.  The  fine  grid  surrounding  this  region  has 
to  be  continued  in  this  refinement  region.  The  numerical  resolution  is  given  by  the  length 
of  the  shortest  element  edge. 


wall  boundaries.  A  similar  behavior  can  be  seen  in  the  test  case  of  Straka  et  al.  in  figure 
3.  These  vortices  are  almost  independent  of  the  numerical  resolution  and  cannot  be  used 
for  comparing  the  accuracy  of  different  simulations. 

5.2  Adaptive  versus  Uniform 

With  the  criterion  introduced  in  the  previous  subsection  we  now  compare  a  simulation 
performed  on  an  adaptively  refined  triangular  mesh  with  a  simulation  on  a  uniform  grid. 
The  result  is  shown  in  figure  6.  The  onset  of  the  Kelvin- Helmholtz  instability  and  even  the 
small  vortices  produced  by  the  instability  are  almost  identical  in  both  simulations.  This 
shows  that  the  locally  refined  mesh  is  able  to  produce  the  same  accuracy  as  a  uniform 
mesh. 

But  there  is  one  important  difference  between  these  two  simulations:  the  simulation 
on  the  locally  refined  mesh  is  four  times  faster  than  the  simulation  using  the  uniform 
mesh  (see  table  1).  Only  5%  of  the  CPU  time  is  used  for  adapting  the  grid  and  only  one 
quarter  of  the  elements  is  used,  i.e.  25%  memory  requirements  compared  to  the  uniform 
grid.  This  demonstrates  very  clearly  the  advantage  of  using  adaptive  mesh  refinement  for 
the  simulation  of  atmospheric  convection. 
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Figure  6:  Comparison  between  a  simulation  using  an  adaptive  mesh  and  a  simulation 
using  a  uniform  mesh  after  1000  seconds.  The  resolution  of  the  uniform  mesh  is  equal  to 
the  finest  resolution  of  the  adaptive  mesh.  The  length  of  the  shortest  element  edges  is 
11m  in  both  cases.  All  parameters  of  the  rising  bubble  are  identical  to  those  in  figure  4. 


5.3  Size  of  Refinement  Region 

As  shown  in  the  previous  subsection  the  simulation  using  adaptive  mesh  refinement  is 
much  faster  than  a  simulation  using  a  uniform  mesh.  Nevertheless  both  simulations  have 
the  same  accuracy.  Can  we  improve  the  efficiency  by  using  a  smaller  refinement  region 
without  losing  accuracy?  To  answer  this  question  we  repeated  the  adaptive  simulation  of 
figure  6  with  three  slightly  different  refinement  regions  (figure  7). 

The  onset  of  Kelvin-Helmholtz  instability  is  significantly  earlier  when  reducing  the  size 
of  the  refinement  region.  This  indicates  that  the  numerical  errors  are  increasing.  In  the 
case  of  this  rising  warm  air  bubble  it  seems  to  be  better  to  use  the  large  refinement  region 
as  in  figure  6.  Regardless,  the  computational  effort  is  not  much  reduced  by  the  smaller 
refinement  regions.  When  developing  new  refinement  criteria  in  future  research  we  will 
have  to  repeat  this  test  to  ensure  that  the  adaptive  mesh  gives  the  same  accuracy  as  the 
uniform  mesh. 

6  SUMMARY  AND  OUTLOOK 

In  this  paper  we  presented  a  numerical  model  for  solving  the  stratified  compressible 
Euler  equations.  It  uses  a  high-order  discontinuous  Galerkin  method  based  on  triangular 
elements11  in  combination  with  a  semi- implicit  time-integrator16.  This  avoids  the  severe 
time-step  restriction  of  explicit  schemes.  To  our  knowledge,  this  is  the  first  time  that  a 
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adaptive 

uniform 

number  of  elements: 

4668 

16384 

grid  time: 

0.5h 

O.Oli 

total  time: 

9.7h 

38. 7h 

Table  1:  Number  of  elements  and  CPU  time  for  the  two  simulations  shown  in  figure  6  with 
semi- implicit  time-integration.  The  time  values  are  the  time  used  for  reaching  t  =  1000s 
of  the  warm  air  bubble  test  case  and  the  number  of  elements  is  an  average  value  over  the 
whole  time.  The  grid  time  denotes  the  CPU  time  that  was  used  for  adapting  the  grid 
after  each  time  step.  Each  simulation  was  done  on  the  same  single  Linux  CPU.  The  finest 
resolution  in  the  adaptive  simulation  covers  21.7%  of  the  whole  domain. 
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Figure  7:  Adaptive  simulation  as  in  figure  6  for  three  different  refinement  regions  after  900 
seconds.  The  refinement  criterion  is  always  given  by  (24),  but  the  number  of  additional 
fine  elements  surrounding  this  region  is  varied.  From  right  to  left:  one  row  of  fine  elements 
is  added  in  each  simulation.  The  largest  refinement  region  shown  on  the  left  panel  is  the 
same  as  in  figure  6. 


17 


Andreas  Miiller,  Jorn  Behrens,  Francis  X.  Giraldo  and  Volkmar  Wirth 


DG  compressible  Euler  model  has  been  combined  with  an  h-adaptive  mesh  refinement 
for  stratified  flows  that  are  of  interest  in  nonhydrostatic  meteorological  applications.  For 
the  h-adaptivity  we  use  the  function  library  AMATOS23  that  uses  a  very  efficient  space 
filling  curve  approach. 

For  testing  our  numerical  model  we  simulated  three  commonly  used  test  cases.  Our 
results  agree  well  with  the  results  from  the  literature. 

For  considering  the  question  whether  the  results  are  affected  by  using  adaptive  mesh 
refinement,  we  developed  a  new  qualitative  error  measure  for  the  simulation  of  rising  warm 
air  bubbles  by  using  artificial  viscosity  with  a  constant  fairly  large  viscosity  parameter 
fi.  We  found  that  Kelvin-Helmholtz  instability  strongly  depends  on  the  accuracy  of  the 
numerical  simulation.  Therefore  the  time  when  the  instability  becomes  visible  is  a  good 
criterion  for  comparing  different  simulations. 

With  this  criterion  we  compared  a  simulation  using  an  adaptive  mesh  with  a  simulation 
using  a  uniform  mesh;  the  finest  resolution  is  the  same  in  both  cases.  Figure  6  shows  al¬ 
most  no  difference  between  the  two  simulations  even  after  the  instability  has  fully  evolved. 
The  simulation  on  the  locally  refined  mesh  is  as  accurate  as  the  simulation  on  the  uniform 
mesh.  But  the  simulation  using  adaptive  mesh  refinement  is  four  times  faster  than  the 
simulation  on  the  uniform  mesh  and  only  5%  of  the  CPU  time  is  used  for  adapting  the 
grid.  This  indicates  that  adaptive  mesh  refinement  should  be  a  big  advantage  for  our 
future  cloud  simulations. 

For  simulating  clouds  we  are  currently  working  on  including  moisture  into  our  model 
and  on  extending  the  model  to  3D. 
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